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We present a new analysis of J /ijj production yields in deuteron-gold collisions at ^/snn ~ 200 
GeV using data taken by the PHENIX experiment in 2003 and previously published in [S.S. Adler 
et al., Phys. Rev. Lett 96, 012304 (2006)]. The high statistics proton-proton J/ip data taken 
in 2005 is used to improve the baseline measurement and thus construct updated cold nuclear 
matter modification factors {RdAu). A suppression of J/ip in cold nuclear matter is observed as one 
goes forward in rapidity (in the deuteron-going direction), corresponding to a region more sensitive 
to initial state low-x gluons in the gold nucleus. The measured nuclear modification factors are 
compared to theoretical calculations of nuclear shadowing to which a J/ip (or precursor) breakup 
cross section is added. Breakup cross sections of ai,reakup ~ 2.8l]^'4 i'^-'^-i't) rnb are obtained by 
fitting these calculations to the data using two different models of nuclear shadowing. These breakup 
cross-section values are consistent within large uncertainties with the 4.2 ± 0.5 mb determined at 
lower collision energies. Projecting this range of cold nuclear matter effects to copper-copper and 
gold-gold collisions reveals that the current constraints are not sufficient to firmly quantify the 
additional hot nuclear matter effect. 

PACS numbers: 25.75.Dw 



I. INTRODUCTION 

Understanding the behavior of QCD matter under dif- 
ferent conditions of temperature and density is the sub- 
ject of intense experimental and theoretical work in nu- 
clear physics. The transition from hadronic matter to 
a quark-gluon plasma at high temperature is expected 
to be achieved in high energy heavy ion collisions. The 
hadronization of partons in vacuum or cold nuclear mat- 
ter (i.e. in a nucleus) is also of keen interest, and repre- 
sents a nonperturbative and dynamic QCD process. The 
formation and interaction of heavy quarkonia (for exam- 
ple J/ip mesons) in vacuum, cold nuclear matter and hot 
nuclear matter present an excellent laboratory for gaining 
insights on these transformations. Recent results from 
the Relativistic Heavy Ion Collider reveal a significant 
suppression of the final J/^j yield in central (small im- 
pact parameter) Au + Au reactions at y^sjvw=200 GeV, 
relative to expectations scaled from p -{- p reactions at 
the same energy [ll, [2[ . A possible source of this suppres- 
sion is the screening of the attractive interaction between 
the quark-antiquark pair in the hot nuclear medium, as 
temperatures are expected to be above the critical tem- 
perature for a quark-gluon plasma transition. Larger 
J/f/) suppression is observed at forward rapidity than at 
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midrapidity, which contradicts models with only color 
screening of quarkonia proportional to the local energy 
density. 

Produced cc pairs must pass through the remaining nu- 
clear material from the incident cold nuclei, in addition 
to surviving any hot medium environment. The so-called 
cold nuclear matter effects 3], including modification of 
initial parton distribution functions (shadowing, gluon 
saturation, anti-shadowing, EMC effect, etc.), initial and 
final state partonic multiple scattering, and related initial 
state parton energy loss need to be accounted for before 
firm conclusions can be drawn about the effect of the hot 
medium thought to be created. In fact, these various cold 
nuclear matter effects are interesting in their own right, 
notably in terms of hadronization time scales, parton en- 
ergy loss in matter, and the various initial state effects 
mentioned above. 

This paper presents a new analysis of the modification 
of J/ip production in deuteron-gold {d -\- Au) collisions 
relative to proton-proton {p-\-p) collisions at ^snn~200 
GeV and the implications for understanding the Au + Au 
and Cu -t- Cu data at the same energy. The PHENIX ex- 
periment has previously published a result using p + p 
and d -\- Au data taken in 2003 'i\. A modest J/tp sup- 
pression was observed at forward rapidity (i.e. in the 
deuteron moving direction), which is a possible indica- 
tion of shadowing of low- a; gluons in the gold nucleus. A 
substantially larger (more than an order of magnitude) 
p -\- p data set was recorded in 2005 with the J/-0 re- 
sults published in 0, and has been used as the baseline 



for recent Au + Au and Cu + Cu nuclear modification 
factors [l|, Q . The same p + p data set is used in the 
analysis presented here to determine the d + Au nuclear 
modifications more accurately and in a fully consistent 
way with those in the Au + Au and Cu + Cu cases. In 
addition, during the two years between the analyses of 
the 2003 and 2005 data sets significant improvements in 
the reconstruction software and signal extraction method 
were achieved, as well as an overall better understanding 
of the detector performance. These improvements have 
been included in this analysis, allowing maximal cancel- 
lation of systematic errors when using the 2005 p+p data 
sample to form the J/tp nuclear modification factor. We 
first describe the updated analysis, then present the new 
nuclear modification factors and their implications. 



II. EXPERIMENT 

The PHENIX apparatus is described in 6] . It consists 
of two sets of spectrometers referred to as the central 
arms, which measure particles emitted at midrapidity 
{\y\ < 0.35), and the muon arms, measuring particles 
emitted at backward and forward rapidity (—2.2 < y < 
-1.2 and 1.2 < y < 2.2). 

At midrapidity, J/^p particles are measured via their 
decay into two electrons. Electrons are identified by 
matching tracks reconstructed with drift chambers (DC) 
and pad chambers (PC) to clusters in the Electromag- 
netic Calorimeters (EMCAL) and hits in the Ring Imag- 
ing Cerenkov Counters (RICH). In d + Au collisions, a 
charged track is identified as an electron candidate by 
requiring at least three matching RICH phototube hits 
within a certain radius with respect to the the center 
defined by the track projection at the RICH, a position 
matching of ±4 standard deviations between the EM- 
CAL cluster and the reconstructed track, and a cut on 
the ratio of energy to momentum. In p -\- p collisions the 
electron identification cuts are the same except that only 
two matching RICH phototube hits are required. 

At forward and backward rapidity, J/ip particles are 
measured via their decay into two muons. Muons are 
identified by matching tracks measured in Cathode Strip 
Chambers (referred to as the Muon Tracker, or MuTR) 
to hits in alternating planes of larocci tubes and steel 
absorbers (referred to as the Muon Identifier, or MuID). 
Each muon arm is located behind a thick copper and iron 
absorber that is meant to absorb most hadrons produced 
during the collisions, so that the measured muons must 
penetrate 8 to 11 interaction lengths of material in total. 

The d + Au data used for this analysis were recorded 
in 2003 using a minimum bias trigger that required hits 
in each of the two beam-beam counters (BBC) located 
at positive and negative rapidity (3 < \ri\ < 3.9), and 
represent integrated luminosities for the different spec- 
trometers ranging from 1.4 nb~^ to 1.7 nb~^ (or equiva- 
lently 2.7 to 3.4 billion interactions). This trigger covers 
88 ± 4% of the total d + Au inelastic cross section of 2260 



mb y\. For the electrons, an additional trigger was used 
that required one hit above threshold in the EMCAL and 
a matching hit in the RICH. For the muons, two addi- 
tional triggers were used at different times during the 
data-taking period. The muon triggers are based on in- 
formation from the MuID, which has five active detector 
layers between the steel absorbers. For the first part of 
the data-taking period, one of the tracks was required to 
reach the fourth MuID plane, while the other was only 
required to reach the second MuID plane. For the latter 
part, the trigger required at least two tracks to reach the 
fourth MuID plane of larocci tubes. 

The BBCs are also used to determine the centrality of 
the d + Au collisions by measuring the energy deposited 
in the counters located at negative rapidity (in the gold- 
going direction). For a given centrality bin, the average 
number of equivalent nucleon-nucleon collisions (iVcoii) is 
derived from this energy using a Glauber calculation i8| 
coupled to a simulation of the BBC. The centrality bins 
used in this analysis and the corresponding number of 
collisions are listed in Table HI To ensure that the cen- 
trality categories are well defined, collisions are required 
to be within ±30 cm of the center of the interaction re- 
gion. 



III. SIGNAL EXTRACTION 

The number of J/ip particles is determined using the 
invariant mass distribution of unlike-sign lepton pairs. 
At midrapidity, the J/V' signal count is obtained via 
counting the number of unlike-sign dielectrons after sub- 
tracting the like-sign pairs in a fixed mass window 2.6 < 
Me+e- < 3.6 GeV/c2 or 2.7 < M^+e- < 3.5 GeV/c^, 
depending on the number of DC hits required for track 
reconstruction. Figure [T] shows the J/ip mass spectrum 
after subtracting the background. The solid black line 
is the sum of the J/ip line shape (dashed curve) and an 
exponential function (dot-dashed curve) describing the 
continuum component determined from the 2005 p + p 
data set Q]. The J/i/; line shape function takes into 
account the momentum resolution of track reconstruc- 
tion, internal radiative effects [9|, and external radiative 
effects evaluated using a GEANT [l^] simulation of the 
PHENIX detector. The number of J/ip particles in d+Au 
collisions is too small to allow a good fit but a compari- 
son between d-l- Au and p+p J/ip line shapes shows good 
agreement. The fraction of J/ip candidates outside of the 
mass window due to the radiative effects is estimated to 
be 7.2% ± 1.0% based on the line shape functions. The 
J/ip signal is also corrected for the dielectron continuum 
yield, which originates primarily from open charm and 
Drell-Yan pairs inside the mass window. The estimated 
contribution is 10% ± 5%, based on the fitting function 
and PYTHIA ^] simulations. Approximately 400 J/ip 
mesons are obtained. 

At backward and forward rapidity an event mixing 
technique is now used to estimate the combinatorial back- 



TABLE I; Characterization of the colhsion centrahty for d + Au colhsions. First hne: centrahty bins used in this analysis; 



Second hne: corresponding number of binary colhsions A'^coii; 
the three rapidity ranges used for this analysis. 



lines 3, 4 and 5: 



3(cent)/f.j/v, for J'/V' mesons emitted 



centrahty 


0-20 % 


20-40 % 


40-60 % 


60-88 % 


0-100 % 


iVcoll 


15.4 ±1.0 


10.6 ±0.7 


7.0 ±0.6 


3.1 ±0.3 


7.6 ±0.3 


c {\y\ < 0.35) 


0.95 ± 0.03 


0.99 ±0.01 


1.03 ±0.01 


1.04 ± 0.03 


0.94 ± 0.02 


c(-2.2 <y < -1.2) 


0.93 ± 0.03 


0.99 ±0.01 


1.04 ±0.01 


1.05 ±0.03 


0.94 ±0.02 


c(1.2 <y < 2.2) 


1.00 ±0.03 


1.00 ±0.01 


1.02 ±0.01 


1.02 ±0.03 


0.94 ± 0.02 



ground, whereas the like-sign pairs were used in the pre- 
vious analysis [J|. The previous method suffered from 
a larger statistical uncertainty for bins where the signal 
over background is poor. A sample mass distribution 
after the subtraction of the mixed event background is 
shown in Figure [21 Approximately 500 and 750 J/ip 
mesons are obtained for backward and forward rapid- 
ity, respectively. The signal counts are determined from 
this subtracted dimuon invariant mass distribution with 
a log-likelihood fit and assuming three different func- 
tional forms and parameters. In all three cases, an ex- 
ponential form is used to account for correlated physical 
background sources (e.g. Drell-Yan or open charm) and 
the possible systematic offset in the normalization of the 
mixed event background. The number of J/ijj particles is 
then estimated either by direct counting of the remain- 
ing number of pairs above the exponential in the mass 
range 2.6 < M^+^- < 3.6 GeV/c^; using a Gaussian 
function with the center fixed to the J/ip mass and the 
width and integrated yield left free; or using two Gaus- 
sian functions for which both the center and widths are 
fixed to the values measured in p + p collisions. The two 
Gaussian functions account for the nonGaussian tails in 
the invariant mass distribution. The normalization of 
the mixed background is varied by a systematic uncer- 
tainty of ±2% prior to its subtraction from the mass dis- 
tribution. This uncertainty is determined by comparing 
different normalization methods. The corresponding sig- 
nal variations are included in the systematic uncertainty. 
Due to the fit procedure described above, for all p+p and 
d + Au cases this normalization uncertainty results in a 
very small systematic uncertainty on the number of mea- 
sured J/ip particles. This entire procedure is identical to 
the one used in [i|, Q . 



IV. INVARIANT YIELD 

The J/ip invariant yield in a given centrahty, trans- 
verse momentum and rapidity bin is: 



Bu <PN, 



J/^ 



N. 



J/V 



^BBC 

"MBicent) 



BBC 



27rpT dpTdy 2TTpTApTAy Aerec^trigN^^f c^^^ 

(1) 

with Bu being the branching ratio for J/^ — > I'^l ; 
-^counts the number of reconstructed J/^p mesons; N^^^ 
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FIG. 1: (color online) Invariant mass spectra in minimum bias 
d + Au reactions for J/ip — > e'^e" at \y\ < 0.35, with the 
functional forms used to extract the number of reconstructed 
J/ip mesons. 



the number of minimum bias events sampled; ^^.psicent) 
the BBC trigger efficiency for minimum bias events in 
a given centrahty category; Apx and Ay the pt and y 
bin widths; A and e^ec the acceptance and reconstruction 
efficiency corrections; etrig the additional J/ip trigger effi- 
ciency and £j^ the BBC efficiency for events containing 
a J/ip. All invariant yields as a function of pt and y in- 
cluding statistical and systematic uncertainties are given 
in Table ni 

The experiment measures the number of J/ijj particles 
per BBC triggered events, which in d + Au collisions rep- 
resent only 88 ± 4% of the total inelastic cross section. 
An additional correction is then applied such that the in- 
variant yield represents 100% of the total inelastic cross 



TABLE 11: The errors quoted in these tables represent 1) the statistical and Type A systematic uncertainties added in quadra- 
ture, and 2) the Type B systematic uncertainties. The Type C errors are shown in the relevant figures. 

Backward Rapidity Results 



Ccntrality (%) 


PT (GcV/c) 


y 




Invariant Yield 








RdA-^ 


0-100 


All 


-2.2, 


-1.2] 


(4.264 ±0.326 ±0.923) 


X 


10" 


-6 


0.90 ±0.08 ±0.19 


0-100 


All 


-2.2, 


-1.7] 


(3.583 ±0.395 ±0.775) 


X 


10" 


-6 


0.95 ±0.12 ±0.20 


0-100 


All 


-1.7, 


-1.2] 


(5.292 ±0.483 ±1.145) 


X 


10" 


-6 


0.90 ±0.09 ±0.19 


0-100 


0-1 


-2.2, 


-1.2] 


(3.040 ±0.460 ±0.658) 


X 


10" 


"^ 


0.69 ±0.11 ±0.15 


0-100 


1-2 


-2.2, 


-1.2] 


(1.782 ±0.201 ±0.386) 


X 


10" 


"^ 


0.84 ±0.10 ±0.18 


0-100 


2-3 


-2.2, 


-1.2] 


(8.141 ± 0.937 ± 1.762) 


X 


10" 


-S 


1.44 ±0.18 ±0.31 


0-100 


3-4 


-2.2, 


-1.2] 


(1.789 ±0.359 ±0.387) 


X 


10" 


-8 


1.21 ±0.26 ±0.26 


0-100 


4-5 


-2.2, 


-1.2] 


(4.016 ± 1.451 ±0.869) 


X 


10" 


-9 


1.14 ±0.43 ±0.24 


0-20 


All 


-2.2, 


-1.2] 


(9.084 ± 0.922 ± 1.925) 


X 


10" 


-6 


0.94 ±0.10 ±0.21 


20-40 


All 


-2.2, 


-1.2] 


(3.676 ±0.642 ±0.770) 


X 


10" 


-6 


0.55 ±0.10 ±0.12 


40-60 


All 


-2.2, 


-1.2] 


(4.013 ±0.583 ±0.842) 


X 


10" 


-6 


0.92 ±0.14 ±0.21 


60-88 


All 


-2.2. 


-1.2] 


(2.062 ±0.312 ±0.436) 


X 


10" 


-6 


1.07 ±0.17 ±0.25 



Mid-rapidity Results 



Ccntrality (%) 


PT (GeV/c) 


V 


Invariant Yield 


RdAu 


0-100 


All 


[-0.35,0.35] 


(6.750 ± 0.540 ± 0.950) x 10"" 


0.85 ±0.07 ±0.15 


0-100 


0-1 


[-0.35,0.35] 


(6.700 ± 0.800 ± 0.940) x 10"^ 


1.05 ±0.14 ±0.21 


0-100 


1-2 


[-0.35,0.35] 


(2.400 ± 0.340 ± 0.340) x 10"^ 


0.74 ±0.11 ±0.15 


0-100 


2-3 


[-0.35,0.35] 


(1.200 ± 0.190 ± 0.170) X 10"^ 


0.96 ±0.17 ±0.19 


0-100 


3-4 


[-0.35,0.35] 


1.37 X 10"** (90% CL) 


0.41 (90% CL) 


0-100 


4-5 


[-0.35,0.35] 


(7.500 ± 3.600 ± 1.100) x lO"** 


1.09 ±0.61 ±0.22 


0-20 


All 


[-0.35,0.35] 


(1.144 ± 0.160 ± 0.160) X 10"^ 


0.71 ±0.10 ±0.12 


20-40 


All 


[-0.35,0.35] 


(7.990 ± 1.290 ± 1.120) x 10"'^ 


0.71 ±0.12 ±0.11 


40-60 


All 


[-0.35,0.35] 


(6.800 ± 1.010 ± 0.950) x lO-"* 


0.93 ±0.14 ±0.14 


60-88 


All 


[-0.35,0.35] 


(3.030 ± 0.500 ± 0.420) x 10"" 


0.94 ±0.16 ±0.14 



Ccntrality (%) 



PT (GeV/c) 



Forward R,apidity Results 

y 



Invariant Yield 



RdA 



0-100 
0-100 
0-100 
0-100 
0-100 
0-100 
0-100 
0-100 
0-20 
20-40 
40-60 
60-88 



All 
All 
All 
0-1 
1-2 
2-3 
3-4 
4-5 
All 
All 
All 
All 



[1.2,2.2] 
[1.2,1.7] 
[1.7,2.2] 
[1.2,2.2] 
[1.2,2.2] 
[1.2,2.2] 
[1.2,2.2] 
[1.2,2.2] 
[1.2,2.2] 
[1.2,2.2] 
[1.2,2.2] 
[1.2,2.2] 



(3.300 
(4.522 
(2.406 
(2.779 
(1.362 
(4.667 
(1.472 
(2.842 
(5.705 
(4.577 
(2.950 
(1.671 



±0.242 
±0.341 
± 0.224 
±0.285 
±0.115 
±0.566 
±0.225 
± 0.756 
±0.501 
± 0.474 
±0.347 
±0.195 



±0.592 
±0.811 
± 0.432 
± 0.498 
± 0.244 
±0.837 
± 0.264 
±0.510 
± 0.987 
± 0.783 
±0.505 
±0.289 



X lO-"* 
X 10"" 
X 10"" 
X 10"'' 
X lO"'^ 
X 10"** 
X 10"** 
X 10"3 
X 10"" 
X 10"" 
X 10"" 
X 10"" 



0.63 
0.68 
0.59 
0.55 
0.60 
0.73 
0.93 
0.84 
0.54 
0.62 
0.62 
0.79 



±0.06 
±0.06 
±0.06 
±0.06 
±0.06 
±0.10 
±0.16 
±0.25 
±0.05 
±0.07 
±0.08 
±0.10 



±0.11 
±0.11 
±0.10 
±0.09 
±0.10 
±0.12 
±0.16 
±0.14 
±0.09 
±0.11 
±0.11 
±0.15 



section (as done in previous PHENIX d + Au analyses). 



The correction factor ratio e 



BBC 



BBC 



MB(cent)/^J/ip 



depends a 



priori on the centrality bin and the rapidity range of the 
measured J/ip particles. The values are given in Table HI 
The same procedure is applied for p+p collisions, so that 
the yields are normalized to the p+p total inelastic cross 
section of 42 mb. 

The acceptance and efficiency corrections are deter- 
mined using a full GEANT simulation [10| of the detector 
with realistic resolutions and detector plane efficiencies 
determined using real data. Compared to the original 
result 4], this simulation benefits from improvements in 
the understanding of the detector alignment, resolution. 



and overall performance. It also includes the improve- 
ments added to the reconstruction software and used for 
the recent p -I- p, Cu -|- Cu and Au + Au analyses [l|,y,la]. 
Although the additional underlying hit occupancies per 
event are modest in p+p and d + Au collisions, they are 
accounted for by embedding the simulated J/ip mesons 
in real data events. The observed differences (4-5%) be- 
tween embedded and nonembedded events are not signif- 
icant given the statistics of the simulations, and therefore 
are included only as a contribution to the systematic un- 
certainty. 

The systematic uncertainties in the J/ip invariant yield 
(Table IIIip are grouped into three categories as in the 
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M„.„ (GeV/c^) 
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FIG. 2: (color online) Invariant mass spectra in minimum 
bias d + Au reactions for (left) J/ijj — > /i^/i~ at —2.2 < y < 
-1.2 and (right) J/tp — > fj,+fi~ at 1.2 < y < 2.2, with the 
functional forms used to extract the number of reconstructed 
J/ip mesons. 



TABLE III: Sources of systematic uncertainties on the J/ip 
invariant yield in d + Au collisions. Columns 2 (3) are the 
average values at mid (forward) rapidity. When two values 
are given, the first (second) is for peripheral (central) colli- 
sions. Uncertainties of type A (type B) are point to point 
uncorrelated (correlated). 



source 


iyj < 0.35 


|y|e [1.2,2.2] 


type 


signal extraction 


6% 


<10% 


A 


acceptance 


8% 


10% 


B 


efficiency 


6% 


8 to 20 % 


B 


run by run variation 


5% 


8% 


B 


input y, pt distributions 


2% 


4% 


B 


embedding 


4% 


5% 


B 



previous analyses: point to point uncorrelated (Type A), 
for which the points can move independently from one 
another; point to point correlated (Type B), for which 
the points can move coherently though not necessarily 
by the same amount; and global uncertainties (Type C), 
for which all points move by the same multiplicative fac- 
tor. Statistical and uncorrelated systematic uncertainties 
(Type A) are summed in quadrature and represented as 
vertical bars. Type B uncertainties are represented with 
boxes. The Type C globally correlated systematic uncer- 
tainties are quoted directly on the Figures. 

Figure [3] shows the invariant J/ip yield as a function 
of transverse momentum for d + Au collisions from this 
new analysis together with the published invariant yield 
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FIG. 3: (color on-line) J/^ invariant yield versus pt in d -I- 
Au collisions and p + p collisions. The three panels are for 
rapidity selections —2.2 < y < —1.2, \y\ < 0.35, and 1.2 < 
y < 2.2 from top to bottom. See text for description of the 
uncertainties and details of the functional fits. 



measured in p+p collisions [2] . From these yields, a (pt^) 
is calculated using the following generic functional form 
to fit the data: 



PrdpT 



A{i + {pT/Bfy 



(2) 



In order to account for finite px binning, the fit func- 
tion is first integrated over each Apx range and the inte- 
gral is compared to data in the corresponding bin. The 
measured (pt^) as well as the associated statistical and 
systematic uncertainties are shown in Table IIVI 

In previous J/ip analyses [ll |2| , it was found that only 
for the high statistics p + p data set (where the measure- 
ment has good precision out to px ~ 8 GeV/c) is the 



8 



TABLE IV: (pt^) calculated from a fit to the data and re- 
stricted to the range < pt < 5 GeV/c. See text for descrip- 
tion of the uncertainties. 



species 



rapidity 



(pt"}[0,5] 



d-l-Au 
d-l-Au 
d-l-Au 



[-2.2,-1.2] 
-0.35, 0.35] 

[1.2,2.2] 



4.3 ±0.3 ±0.4 
3.9 ±0.3 ±0.3 
4.0 ±0.2 ±0.4 



p±p 
P±P 
P±P 



[-2.2,-1.2] 
-0.35, 0.35] 

[1.2,2.2] 



3.4±0.1±0.1 
4.1 ±0.2 ±0.1 
3.4 ±0.1 ±0.1 



functional form of the px spectrtim well constrained. In 
the Au + Au case, the functional form is not well con- 
strained and leads to a very large systematic uncertainty 
on the (pt^) if integrated from to oo. The integral was 
therefore limited to px < 5 GeV/c, where it is best con- 
strained by the data. The d + Au data set suffers from 
the same statistical limitations and the same truncation 
to pt < 5 GeV/c is applied. Finally, this constraint is 
also applied to the p+p case to make a direct comparison 
possible. 

Two uncertainties are quoted in Table IIVI The first 
corresponds to the statistical and point-to-point uncorre- 
lated systematic uncertainties (Type A) on the measured 
yields. It is obtained directly from the fit using the sec- 
ond derivatives of the x^ surface at the minimum. The 
second corresponds to the point-to-point correlated sys- 
tematic uncertainties (Type B). The contribution from 
the Type B uncertainty is estimated independently by 
coherently moving the measured points within the one 
standard deviation limit given by these uncertainties, al- 
lowing them to be either correlated or anti-correlated, 
and then re-doing the fit in all cases. The largest differ- 
ence observed in the values obtained by the fit is used as 
an upper limit to the 1-sigma point-to-point correlated 
uncertainties on the (pt^)- 

In the previous publication [Ij , values for the fully inte- 
grated (pt^) in p+p and d+Au are quoted. A significant 
systematic uncertainty originating from not knowing the 
functional form to best describe the data was found since 
then that was not included in the uncertainty quoted in 
the paper. In addition, the new analysis revealed a bias 
in the previous result that increased the signal, particu- 
larly in the lowest px bin. This bias is now corrected by 
using the mixed event background subtraction technique 
described above together with the modified log-likelihood 
fit over a more appropriate range, corresponding to the 
region where the physical background can accurately be 
described by a single exponential function. Finally, no 
separate treatment of the point-to-point correlated sys- 
tematic uncertainties was performed at that time, since 
it was assumed that it would move all points in the same 
direction (positive correlation) and thus have no impact 
on the measured (pt^)- 

The data, within uncertainties, includes the possib- 



lity of a modest broadening of the transverse momentum 
distribution relative to p + p collisions. This is often 
attributed to initial and final state multiple scattering, 
sometimes referred to as the "Cronin effect." However, 
in calculating the A(pT^) = {pT^)dAu — (PT^)pp one finds 
this effect needs reduced uncertainties from future larger 
data sets to make any firm conclusions. 

Figure [4] shows the J/ip invariant yield, integrated over 
all pt, as a function of rapidity for d + Au collisions. 
Shown are the results of the new analysis presented in 
this paper, as well as the previously published results [jj 
using the same data set. Overall the agreement of the 
two analysis results is good. The two sets of points dif- 
fer in the reconstruction software, analysis cuts, and sig- 
nal extraction technique. Thus many of the systematic 
uncertainties are different, and even the statistical un- 
certainties are not identical due to the different analysis 
cuts and the use of event mixing to estimate the combi- 
natorial background in the new analysis, as opposed to 
the like-sign mass distribution used in [4]. 

Figure [5] shows the J ftp invariant yield for p + p colli- 
sions, from both the published high statistics result from 
Run-5 [3| , as well as the lower statistics result from Run-3 
as published in [^. In both cases the points are in good 
agreement within the systematic uncertainty bands. A 
new analysis of the Run-3 p+p lower statistics data set 
was also performed using the same technique and anal- 
ysis cuts as for d + Au collisions. It also shows good 
agreement with these two sets of measurements, albeit 
with larger statistical uncertainties. 
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FIG. 4: (color online) J/tp invariant yield as a function of 
rapidity for d + Au collisions. Shown are the new analysis re- 
sults from this paper, in addition to the originally published 
results [J] using the same data. The global systematic uncer- 
tainty quoted is for the new analysis. 
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FIG. 5: (color online) J/t/j invariant yield as a function of 
rapidity for p + p collisions. Shown are the high statistics 
results from 2005 p + p PHENIX data taking period 2] , and 
the originally published results [j using the 2003 p + p data 
set. The global systematic uncertainty quoted is for the new 
analysis. 



In Figures S] and [51 the highest rapidity point is not 
located exactly at the same rapidity position between 
the original and the new analysis. This is due to the 
fact that the positive rapidity muon arm has a shghtly 
larger rapidity coverage than the negative rapidity arm. 
This property was used in the 2003 analysis to include 
additional J/ip mesons at forward rapidity in order to 
probe a slightly lower region of x. It was found, however, 
that there were very few counts in this region and that 
the asymmetric rapidity range created additional difhcul- 
ties when comparing the results measured at forward and 
backward rapidity (in case of symmetric collisions) and 
when comparing the results obtained in p+p collisions to 
Cu + Cu or Au + Au collisions, for which this extra ra- 
pidity coverage was not available (due to high occupancy 
limitations at forward rapidity). As a consequence, it was 
decided for the later analyses to forgo the extra few J/ip 
counts at very forward rapidity and use the same width 
rapidity bins at both positive and negative rapidity. 



V. NUCLEAR MODIFICATION FACTOR 

The J/f/) nuclear modification factor in a given central- 
ity and rapidity bin is: 



RdAu 



1 rfJV^;/7dy 
(A^cou) dNP+P/dy 



(3) 



with dN'^^^ /dy being the J/tl) invariant yield measured 
in d + Au collisions; dNj^,/dy the J/ip invariant yield 
measured in p+p collisions for the same rapidity bin and 
(Ncoii) the average number of binary collisions in the cen- 
trality bin under consideration, as listed in Table H] All 
i?(iAu values as a function of p^, y and centrality includ- 
ing statistical and systematic uncertainties are given in 
the Appendix Table |TT1 

Figure \E\ shows the nuclear modification factor RciAm 
calculated using the d + Au new analysis presented in 
this paper for the numerator and the 2005 p + p data for 
the denominator. In contrast to the previous analysis [J], 
where the p + p results were symmetrized around y = 
before calculating ii^Au to compensate for lower J/ip 
statistics in the 2003 p+p data set, in this case the i?dAu 
values are calculated independently at each rapidity. 

The understanding of the detector performance in 
terms of alignment, resolution and efficiency has signif- 
icantly improved between this analysis and previously 
published PHENIX d + Au results Q- This resulted in 
changes in the reconstruction software, analysis cuts, sig- 
nal extraction technique and handling of both the phys- 
ical and combinatorial background in the dilepton in- 
variant mass distribution. Simultaneously, the system- 
atic uncertainties associated with the measurement have 
also been re-evaluated in a way consistent with what was 
learned for the p + p, Cu -I- Cu, and Au + Au analyses. 
The new uncertainties are in general larger, although 
some of them cancel with their p + p counterpart when 
forming i^^Au. This approximately counterbalances the 
reduction of the statistical uncertainty achieved by us- 
ing the 2005 p + p data set as a reference. Additionally, 
the J/ijj production cross sections in p + p collisions mea- 
sured in 2005 [4I are compatible within uncertainties, but 
higher than the values used in [4l] (based on the 2003 p+p 
data set) by about 13 %. As a consequence, the new nu- 
clear modification factors are systematically lower than 
the ones previously published by about 5 to 20% for most 
points, depending on the pt, y or centrality bin that is 
considered. 

Within uncertainties, the nuclear modification factors 
are consistent with -R^Au = 1.0 at negative and midra- 
pidities, and are significantly lower than 1.0 at forward 
rapidity only, that is, in the deuteron-going direction. 
This trend is similar to that shown in Figure 1 of Q, al- 
though the new values are systematically smaller for all 
rapidity bins. 

Figure [7] shows the J/^p nuclear modification factor 
in d -|- Au collisions as a function of the number of bi- 
nary collisions for three rapidity ranges and four central- 
ity classes. Only at forward rapidity is there statistically 
significant suppression. 



VI. DISCUSSION 

As stated in the Introduction, the d + Au data is 
interesting both to fundamentally understand issues of 
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FIG. 6: J/ijj nuclear modification factor -RjAu as a function 
of rapidity. 



quarkonia and cold nuclear matter and also to separate 
these effects from hot nuclear matter effects in heavy ion 
collisions. In order to address both issues, we compare 
the experimental data with two different models includ- 
ing both modification of the initial parton distribution 
functions (PDF) and a free parameter to account for the 
breakup of correlated cc pairs that might have otherwise 
formed J/-0 mesons. Note that often in the literature, 
this breakup process in cold nuclear matter is referred to 
as an absorption cross section of the J /ip particles on the 
nucleons in the nucleus. Here we avoid this nomenclature 
both because the object that is "absorbed" is generally 
not a fully- formed J/?/; but rather a cc pair, and because 
the actual process is more a breakup of this pair, rather 
than the absorption of it. 

Shown in Figure [8] is the nuclear modification factor 
-R^Au as a function of rapidity in comparison to theo- 
retical calculations 3] that include either EKS ^] or 
NDSG ^\^ shadowing models for the nuclear PDFs. In 
each case an additional suppression associated with a 
(^breakup IS also included. Note that there is no ah ini- 
tio calculation of this cross section, and while one might 
expect a similar value to results at lower energy ^14], it 
need not be identical. 

Taking full account of the statistical and systematic 
uncertainties on the experimental data, the breakup cross 
section is determined under certain assumptions. We 
have followed the statistical procedure detailed in [15j . 
If we assume that the EKS modified nuclear PDFs are 
exactly correct, and that the only additional suppression 
is accounted for by Ubreakup, then the data constrains 
^breakup = '^■&-{4 mb with the uncertainties as one stan- 
dard deviation. Similarly, if we assume the NDSG modi- 
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FIG. 7; J/V" nuclear modification factor RdAn as a function 
of A'coii for three rapidity ranges. 



fied nuclear PDFs, then we obtain abreakup = 2.2^i^ mb. 
These breakup cross-section values are consistent (within 
the large uncertainties) with the 4.2 ±0.5 mb determined 
at lower energies at the CERN-SPS 14]. The extracted 
breakup cross section at lower energies assumes no con- 
tribution from the modification of nuclear PDFs. At the 
lower energies, J/^ production is sensitive to higher- a; 
partons in the anti-shadowing regime where the modifi- 
cations are expected to be smaller and in the opposite 
direction [iq . 

The modified nuclear PDFs from EKS and NDSG are 
constrained from other experimental measurements such 
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FIG. 8: (color online) RdAn data compared to various the- 
oretical curves for diflerent abreakup values. Also, shown as 
a band are the range of abreakup found to be consistent with 
the data within one standard deviation. The top panel is a 
comparison for EKS shadowing [I4], while the bottom panel 
is for NDSG shadowing [11]. 



FIG. 9; (color online) RdAu data as a function of A^coii 
for three different rapidity ranges. Overlayed are theoreti- 
cal curves representing the best fit ^breakup values as deter- 
mined in each rapidity range separately, utilizing EKS and 
NDSG nuclear PDFs and a simple geometric dependence. 
Also, shown as bands are the range of abreakup found to be 
consistent with the data within one standard deviation. 



as deep inelastic scattering from various nuclear targets 
and the resulting F2{A) structure functions. A geometric 
parametrization of these PDFs based on the path of the 
parton through the nucleus is described in 13| and Q. 
One can test this geometric dependence by comparison 
with the d-|- Au nuclear modification factors as a function 
of iVcou. Using this geometric dependence, the most prob- 
able (Jbreakup IS Calculated independently in three rapidity 



ranges (see Table |V]) . The corresponding nuclear modifi- 
cation values and their one standard deviation bands are 
shown as a function of iVcou in Figure [51 The two calcu- 
lations with EKS and NDSG nuclear PDFs yield almost 
identical bands since the same geometric dependence is 
used in both cases. However, each band represents a dif- 
ferent balance of modification due to the nuclear PDF 
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TABLE V: Most probable values and one standard deviations 
of (Jbreakup assuming two different shadowing models, from a 
fit to minimum bias RdAu points as a function of rapidity 
(Figure IS]), and fits to -R^Au as a function of A'^coii in three 
separate rapidity bins (Figure [Oj. 
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and the breakup cross section. 

For both the EKS and the NDSG PDFs, the values of 
Cbreakup extracted from the overall rapidity dependence 
of i?dAu and from the A^coii dependence of RdAu within 
the different rapidity ranges are consistent within the 
large systematic uncertainties. It should be noted that 
though the 1—p- value for the best fit is poor at backward 
rapidity (as can be seen in Fig. [9]), there is still a well de- 
fined maximum in the likelihood function for atreakup- 
A future higher precision Rd-Au measurement as a func- 
tion of centrality will be crucial to constrain the exact 
geometric dependence. 

One can also utilize this model to do a consistent cal- 
culation of the contribution from cold nuclear matter ef- 
fects that should be present in Cu -I- Cu and Au + Au 
collisions. These contributions, using the best fit value 
of (Jbreakup and their one standard deviation values ex- 
tracted from the data in Figure [5] for each of the two 
shadowing models, are shown in Figures [TOl and [TTl In 
the Cu + Cu case, J/tfj production is not suppressed be- 
yond cold nuclear matter effects at midrapidity or at for- 
ward rapidity, within the limits of the large error bands, 
and the midrapidity data in the Au + Au case is similarly 
inconclusive. However, there is a significant suppression 
in the data at forward rapidity, beyond the uncertain- 
ties in both the data and the projection. It should be 
noted that the uncertainty band at forward and midra- 
pidities are entirely correlated, as they reflect only the 
uncertainty in the (Jbreakup- There is no systematic un- 
certainty included for the choice of modified nuclear PDF 
model, which is the only way to change the relative sup- 
pression between forward and midrapidities within the 
context of this calculation. The more data-driven cal- 
culation described later in this section, however, is per- 
formed independently at different rapidities and does not 
suffer the same stipulation. 

It should also be noted that the theoretical calcula- 
tions yield Raa as a continuous function of the number of 
participants, whereas the data points are at discrete val- 
ues representing a convolution of the modification factor 
with the Npart distribution within a particular centrality 
category. A Glauber simulation combined with a Monte 
Carlo of the PHENIX experimental trigger and central- 
ity selection is utilized to convert the continuous the- 
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FIG. 10; (color online) Raa for Cu-|-Cu 0] collisions com- 
pared to a band of theoretical curves for the (Jbreakup values 
found to be consistent with the d + Au data as shown in Fig- 
ure [8] The top figure includes both EKS shadowing [12] and 
NDSG shadowing [l^ at midrapidity. The bottom figure is 
the same at forward rapidity. 



ory predictions into discrete predictions in the simulated 
PHENIX centrality categories. Thus, the results shown 
in the figures are in fact predictions for the matched event 
selection categories of the experimental data points. 

In order to explore the cold nuclear matter constraints 
further, an alternative data-driven method proposed 
in [18] is used. This approach assumes that there is a 
single modification factor parameterizing all cold nuclear 
matter effects that is a simple function of the radial posi- 
tion in the nucleus. This computation has the advantage 
of not having to assume a specific shadowing scheme and 
a specific breakup cross section, but instead relies only on 
the measured impact parameter dependence. It assumes 
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FIG. 11: (color online) Raa for Au+Au [J] collisions com- 
pared to a band of theoretical curves for the abreakup values 
found to be consistent with the d + Au data as shown in Fig- 
ure [S] The top figure includes both EKS shadowing 12] and 
NDSG shadowing 13] at midrapidity. The bottom figure is 
the same at forward rapidity. 



that the cold nuclear matter effects suffered by a J/ip 
in a Au -I- Au collision at a given rapidity are the prod- 
uct of the modifications measured in d -f- Au collisions 
at the same rapidity and the modifications measured at 
the opposite rapidity (or equivalently in a Au -I- d colli- 
sion). This assumption holds for the two effects consid- 
ered so far, namely shadowing and subsequent breakup. 
It also assumes that the same parton distributions are 
sampled by the J/ip particles observed in the (wide) ra- 
pidity range in Au -I- Au and d-|- Au collisions. Note that 
since this model implicitly includes any possible modi- 
fied nuclear PDFs, the modification factors may have an 
x-dependence that is accounted for by considering the 



backward, mid and forward rapidity d + Au data. The 
different rapidity regions are sensitive to the initial-state 
partons in the gold nucleus in three broad ranges of x, 
corresponding to x « 0.002-0.01, 0.01-0.05, and 0.05-0.2, 
as determined from PYTHIA. 

A Glauber Monte Carlo and a simulation of the BBC 
detector used for centrality determination and triggering 
are done. The resulting four centrality categories (0-20%, 
20-40%, 40-60%, 60-88%) in d -f- Au collisions are char- 
acterized by a distribution in the number of binary colli- 
sions, as shown in the top panel of Figure ll2l In addition, 
the distribution of radial positions r in the Au nucleus of 
binary collisions is calculated and shown in the bottom 
panel of Figure [T^l 

The procedure is to use the forward, mid, and back- 
ward rapidity centrality-dependent i?dAu to constrain 
the modification factor 5ft(r) for three broad regions of 
initial parton x (JUiow, ^mid, ^high respectively). Then 
one can use these parameterizations to project the cold 
nuclear matter effect in the Au + Au case. The current 
d + Au data are insufficient to constrain the functional 
form of di{r). As a simplifying case, 5R(r) is assumed to 
be linear in r and to be fixed at 5R(r > 8 fm) = 1.0 at 
the edge of the gold nucleus. Thus, the only free pa- 
rameter is the slope (or equivalently the magnitude of 
the modification factor at r = 0). Other functions were 
tried and essentially differ by their extrapolation to lower 
and higher radial positions, since the data are not precise 
enough to constrain the shape. This has a particularly 
strong impact on the most peripheral collisions for which 
our assumption that ^{r > 8 fm) = 1.0 adds a significant 
constraint to the shape. 

For all possible slope parameters, consistency with the 
experimental data is checked using the procedure de- 
tailed in d5|, which utilizes the full statistical and sys- 
tematic uncertainties. The range of parameters within 
one standard deviation of the uncertainties is deter- 
mined separately for backward, mid, and forward ra- 
pidity. Using this range of parameters, the cold nuclear 
matter suppression expected in Au -I- Au collisions is pro- 
jected as a function of collision centrality and for mid 
and forward rapidity. Note that the forward rapidity 
Au + Au J/ip production is sensitive to the low-x par- 
tons in one gold nucleus and the high-x partons in the 
other gold nucleus. Thus, in the Monte Carlo, for ev- 
ery binary collision at ri and r2 (the radii with respect 
to the center of each nucleus) the expected modification 
is ^lowifi) X ^highir2)- The midrapidity Au + Au J/t/j 
production is predominantly sensitive to the mid- a; par- 
tons from both gold nuclei and therefore the expected 
modification is ^mid{ri) x ^mid{r2)- The total modifica- 
tion expected is calculated by taking the average over 
all correlated ri and r2 positions for binary collisions 
within overall Au + Au collisions in each Au + Au cen- 
trality class. 

The results of these calculations matched to the exper- 
imentally measured Au + Au centrality bins are shown in 
Figure [131 It is notable that the midrapidity cold nuclear 
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matter extrapolation agrees within the uncertainty of the 
experimental data at midrapidity. Thus, it is not possible 
within the current constraints to determine the potential 
extent of hot nuclear matter effects. This conclusion is 
qualitatively similar to that reached from the previous 
model calculations as shown in Figure [TlJ However, at 
forward rapidity, this method projects a somewhat larger 
range of possible cold nuclear matter effects than the pre- 
vious models. 

Neither the predictions of cold nuclear matter effects 
in heavy ion collisions based on fitting of the d + Au data 
with theoretical curves (Figures [TUland[TT|l. nor those ob- 
tained directly from the d + Au data points (Figure [T5|) 
are well enough constrained to permit quantitative con- 
clusions about additional hot nuclear matter effects. 
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FIG. 12: Results from a Glauber Model Monte Carlo includ- 
ing simulation of the d+Au centrality selection and triggering 
based on the PHENIX BBC. The top panel shows the distri- 
bution of the number of binary collisions for events in each 
of the four centrality classes 0-20%, 20-40%, 40-60%, 60-88%. 
The distribution for radial impact points in the gold nucleus 
of binary collisions is shown in the lower panel. 
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FIG. 13: Predictions of the data driven method [I3] con- 
strained by the -R^Au as a function of collision centrality for 
the Au-|-Au Raa for midrapidity (top) and at forward rapid- 
ity (bottom). 



VII. CONCLUSIONS 

A new analysis of J/tp production in d-|-Au collisions at 
Y^snn = 200 GeV has been presented using the PHENIX 
2003 d-|-Au data set. Cuts and analysis techniques that 
are consistent with the previously published results for 
p + PjCu -I- Cu and Au -t- Au collisions at the same en- 
ergy [l|, 0, [3 are used. The new analysis also benefits 
from the significantly larger p-t-p data set from Run-5. 

A statistical comparison of these new results to theo- 
retical calculations has been performed with a detailed 
handling of the experimental uncertainties to estimate a 
J/'0 (or precursor) breakup cross section in cold nuclear 
matter on top of models for the modifications of the par- 
ton distribution functions in the nucleus. Using EKS 
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(NDSG) shadowing, a breakup cross section of 2.8ti;4 
(2.2tJ'5) rnb is obtained. These breakup cross-section 
values are consistent within large uncertainties with the 
4.2 ± 0.5 mb determined at lower energies at the CERN- 
SPS [U]. The measured values are then used to predict 
the expected cold nuclear matter effects on J /ip produc- 
tion in Cu -|- Cu and Au -I- Au collisions, and these are 
compared to the measured nuclear modification factors 
for those systems. These predictions are found to be 
similar to those from a less model-dependent and more 
data-driven method based on the variation of the nu- 
clear modification factor measured in d -I- Au collisions 
as a function of both rapidity and centrality [18|. It is 
notable that the latter method yields a somewhat larger 
possible suppression in the forward rapidity case. In all 
cases the large error bars associated with the extrapo- 
lation prevent making firm quantitative statements on 
any additional J /ip suppression in Au -I- Au collisions be- 
yond that expected from cold nuclear matter effects. A 
d + Au data set with much improved statistical precision 
is needed to both reduce the statistical uncertainties and 
permit better control over the systematic uncertainties. 
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